function B_=B_coil_helmholtz_rho(z_,rho_,I_,a_,d_)
   B_=B_coil_rho(z_+d_,rho_,I_,a_)+B_coil_rho(z_-d_,rho_,I_,a_);
end
